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The tunneling process in a many-body system is a phenomenon which 
hes at the very heart of quantum mechanics. It appears in nature in the 
form of a-decay, fusion and fission in nuclear physics, photoassociation and 
photodissociation in biology and chemistry. A detailed theoretical description 
of the decay process in these systems is a very cumbersome problem, either 
because of very complicated or even unknown interparticle interactions or due 
to a large number of constitutent particles. In this work, we theoretically 
study the phenomenon of quantum many-body tunneling in a more transparent 
and controllable physical system, in an ultracold atomic gas. We analyze a 
full, numerically exact many-body solution of the Schrodinger equation of a 
one-dimensional system with repulsive interactions tunneling to open space. 
We show how the emitted particles dissociate or fragment from the trapped and 
coherent source of bosons: the overall many-particle decay process is a quantum 
interference of single-particle tunneling processes emerging from sources with 
different particle numbers taking place simultaneously. The close relation to 
atom lasers and ionization processes allows us to unveil the great relevance 
of many-body correlations between the emitted and trapped fractions of the 
wavefunction in the respective processes. 

The tunneling process has been a matter of discussion [IH3] since the advent of quantum 
mechanics. In principle, it takes place in all systems whose potential exhibits classically for- 
bidden but energetically allowed regions. See, for example, the overview ref. ^ and Fig. |l| 
When the potential is unbound in one direction, the quantum nature of the systems allows 
them to overcome potential barriers for which they classically would not have sufficient en- 
ergy and, as a result, a fraction of the many-particle system is emitted to open space. For 
example, in fusion, fission, photoassociation and photodissociation processes, the energet- 
ics or life times are of primary interest [5H9]. The physical analysis was made under the 
assumption that the correlation between decay products (i.e., between the remaining and 
emitted fractions of particles) can be neglected. However, it has to be stressed ffist that, at 
any finite decay time, the remaining and emitted particles still constitute one total many- 
body wavefunction and, therefore, can be correlated. Second, in contrast to the tunneling 
of an isolated single particle into open space, which has been amply studied and understood 
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[1], nearly nothing is known about the tunnehng of a many-body system. In the present 
study, we demonstrate that the fundamental many-body aspects of quantum tunneling can 
be studied by monitoring the correlations and coherence in ultracold atomic gases. For this 
purpose, the initial state of an ultracold atomic gas of bosons is prepared coherently in a 
parabolic trapping potential which is subsequently transformed to an open shape allowing 
for tunneling (see the upper panel of Fig. [T]). In this tunneling system the correlation be- 
tween the remaining and emitted particles can be monitored by measuring deviations from 
the initial coherence of the wavefunction. The close relation to atom lasers and ionization 
processes allows us to predict coherence properties of atom lasers and to propose the study 
of ionization processes with tunneling ultracold bosons. 

In the past decade, Bose-Einstein condensates (BECs) [lOl [H] have become a toolbox to 
study theoretical predictions and phenomena experimentally with a high level of precision 
and control [12j. BECs are experimentally tunable, the interparticle interactions [13], trap 
potentials [H], number of particles [IS], their statistics [TB] and even dimensionality [T71 - 
[TU] are under experimental control. From the theoretical point of view, the evolution of 
an ultracold atomic cloud is governed by the time-dependent many-particle Schrodinger 
equation (TDSE) [20] with a known Hamiltonian (details in the Methods section). To study 
the decay scenario in Fig. [T| we solve the TDSE numerically exactly for long propagation 
times. For this purpose, we use the multiconfigurational time-dependent Hartree method 
for bosons (MCTDHB) (for details and literature see the Methods section). Our protocol 
to study the tunneling process of an initially parabolically trapped system into open space 
is schematically depicted in the upper panel of Fig. [T} We restrict our study here to the 
one-dimensional case, which can be achieved experimentally by adjusting the transverse 
confinement appropriately. As a first step, we consider = 2, 4 or 101 weakly repulsive 
^''Rb atoms in the groundstate of a parabolic trap (see the Methods section for a detailed 
description of the considered experimental parameters). 

The initial one-particle density and trap profile are depicted in the upper panel of Fig. [T} 
Next, the potential is abruptly switched to the open form V{x,t > 0) indicated in this 
figure, allowing the many-boson system to escape the trap by tunneling through the barrier 
formed. Eventually, all the bosons escape by tunneling through the barrier, because the 
potential supports no bound states. 

The initial system, i.e., the source of the emitted bosons, is an almost totally coherent 
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state [T] . The final state decays entirely to open space to the right of the barrier, where the 
bosons populate many many-body states, related to Lieb-Liniger states [22H21], which are 
generally not coherent. It is instructive to ask the following guiding questions: what happens 
in between these two extremes of complete coherence and complete incoherence? And how 
does the correlation (coherence) between the emitted particles and the source evolve? By 
finding the answers to these questions we will gain a deeper theoretical understanding of 
many-body tunneling, which is of high relevance for future technologies and applied sciences. 
In particular, this knowledge will allow us to determine whether the studied ultracold atomic 
clouds qualify as candidates for atomic lasers [TT], or as a toolbox for the study of 

ionization or decay processes [SHE]. 

Since the exact many-body wavefunctions are available at any time in our numerical 
treatment, we can quantify and monitor the evolution of the coherence and correlations of the 
whole system as well as between the constituting parts of the evolving wave-packets. In the 
further analysis we use the one-particle density in real p{x, t) and momentum p{k, t) spaces, 
their natural occupations pf^{t) and correlation functions g'^^\k'\k] t), see PSVISU] . We defer 
the details on these quantities to the Methods section. To study the correlation between the 
source and emitted bosons we decompose the one-dimensional space into the internal "IN" 
and external "OUT" regions with respect to the top of the barrier, as illustrated by the 
red line and arrows in the lower panel of Fig. [T] This decomposition of the one-dimensional 
Hilbert space into subspaces allows us to quantify the tunneling process by measuring the 
amount P^ot pi't) of particles remaining in the internal region in real space as a function of 
time. In Fig. [2] we depict the corresponding quantities for N = 2 and = 101 by the green 
dotted curve. A first main observation is that the tunneling of bosonic systems to open 
space resembles an exponential decay process. 

The key features of the dynamics of quantum mechanical systems manifest themselves 
very often in characteristic momenta. Therefore, it is worthwhile to compute and compare 
evolutions of the momentum distributions p{k,t) of our interacting bosonic systems. Fig. [s] 
depicts p{k, t) for A^ = 2, 4, 101 bosons. At t = all the initial real space densities have 
Gaussian-shaped profiles resting in the internal region (see upper panel of Fig. [T]). Therefore, 
their distributions in momentum space are also Gaussian-shaped and centered around /c = 0. 
With time the bosons start to tunnel out of the trap. This manifests itself in the appearance 
of a pronounced peak structure on top of the Gaussian-shaped background, see central upper 
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panel of Fig. |3] The peak structure is very narrow - similar to a laser or an ionization process, 
the bosons seem to be emitted with a very well defined momentum. For longer propagation 
times a larger fraction of bosons is emitted and more intensity is transferred to the peak 
structure from the Gaussian background. Thus, we can relate the growing peak structures 
in the momentum distributions to the emitted bosons and the Gaussian background to 
the bosons in the source. We decompose each momentum distribution into a Gaussian 
background and a peak structure to check the above relation (see Methods for details). The 
integrals over the Gaussian momentum background, Pnot,pi't), are depicted in Fig. |2] as a 
function of time. The close similarity of the Pnot p{^)i characterizing the amount of particles 
remaining in the internal region in real space, and Pnot p{^) confirms our association of the 
Gaussian-shaped background with source bosons and the peak structure in the momentum 
distribution with emitted ones (cf. top panel of Fig. [s]). 

Now we are equipped to look into the mechanism of the many-body physics in the tun- 
neling process with a simplistic model. When the first boson is emitted to open space one 
can estimate its total energy as an energy difference, — E'^~^, of source systems made 
of and — 1 particles. This corresponds to assuming the negligibility of the coupling of 
the trapped and emitted bosons: 

EouT = E'^ — E^ ^ = fii. 

This is simply the chemical potential of the parabolically trapped source system. For non- 
interacting particles its value would be independent of the number of bosons inside the well. 
Since the boson tunnels to open space where the trapping potential is zero, it can be con- 
sidered for the time being as a free particle with all its available energy, /ii, converted to 
kinetic energy Ep^rp = We have to expect the first emitted boson to have the momentum 
ki = yJ2m^i. The value of the estimated momentum agrees excellently with the position of 
the peak in the computed exact momentum distributions, see the arrows marked ki in Fig. [3] 
i)-iv). This agreement allows us to interpret the peak structures in p{k,t) as the momenta 
of the emitted bosons. As a striking feature, also other peaks with smaller k appear in 
these spectra at later tunneling times [see Fig. [3|), ii) and iv)]. We can use an analogous 
argumentation and associate the second peak with the emission of the second boson and 
its chemical potential /i2. Its kinetic energy and the corresponding momentum, k2, can be 
estimated as the energy difference fj,2 = E^~'^ — E'^~^ between the source subsystems made 
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of — 1 and N — 2 bosons, under the assumption of zero interaction between the emitted 
particles and the source. The correctness of the apphed logic can be verified from Fig. |3]i) 
and ii), where for the N = 2 {N = 4) particles the positions of the estimated momenta ^2 
(^3, ^4) of the second (third, fourth) emitted boson fit well with the position of the second 
(third, fourth) peak in the computed spectra. The momenta and details on the calculation 
are collected in the Supplementary Information. The momentum spectrum for = 101 
bosons shows a similar behavior - the multi-peak structures gradually develop with time 
starting from a single-peak to two-peaks and so on, see Fig. [sjiv). However, from this figure 
we see that the positions of the peaks' maxima, and with them the momenta of the emitted 
bosons change with time. On the one hand we see that the considered tunneling bosonic 
systems can not be utilized as an atomic laser: the intially coherent bosonic source emits 
particles with different, weakly time-dependent momenta. In optics such a source would be 
called polychromatic. On the other hand we can associate the peaks with different channels 
of an ionization process and their time-dependency with the channels' coupling. Thus, we 
conclude that it is possible to model and investigate ionization processes with tunneling 
ultracold bosonic systems. 

Let us now investigate the coherence of the tunneling process itself. In the above analysis 
of the momentum spectra we relied on the exact numerical solutions of the TDSE for = 
2,4,101 bosons. We remind the reader that in the context of ultracold atoms the Gross- 
Pitaevskii (GP) theory is a popular and widely used mean-field approximation describing 
systems under the assumption that they stay fully coherent for all times. In our case the 
GP approximation assumes that the ultracold atomic cloud coherently emits the bosons to 
open space and keeps the source and emitted bosons coherent all the time. To learn about 
the coherence properties of the ongoing dynamics it is instructive to compare exact many- 
body solutions of the TDSE with the idealized GP results, see Fig. [sjii). The strengths of 
the inter-boson repulsion have deliberately been chosen such that the GP gives identical 
dynamics for all studied. It is clearly seen that for short initial propagation times the 
dynamics is indeed coherent. The respective momentum spectra obtained at the many-body 
and GP levels are very similar, see Fig. [3j)-iv) for p{k,ti = 100). At longer propagation 
times (if: > ti), however, the spectra become considerably different. This means that with 
time, the process of emission of bosons becomes less coherent. 

Next, to quantify the coherence and correlations between the source and emitted bosons 
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we compute and plot the momentum correlation functions \g^^\k', k\t)\'^ in the left panel of 
Fig. [Ijfor = 101 (for N = 2,A they look almost the same). Let us stress here that the 
proper correlation properties cannot be accounted for by approximate methods. For example 
the GP solution of the problem gives = 1, i.e., full coherence for all times. For the 

exact solution we also obtain that at t = the system is fully coherent, and thus \g^^\k'\k; t = 
0)p = 1. Hence, the top left panel of Fig. |4]is also a plot for the GP time-evolution. However, 
during the tunneling process the many-body evolution of the system becomes incoherent, 
i.e., 0. The coherence is lost only in the momentum-space domain where the 

momentum distributions are peaked, the fc-region associated with the emitted bosons (see 
left panel of Fig. |4]). In the remainder of A;-space the wavefunction stays coherent for all 
times. We conclude that the trapped bosons within the source remain coherent. The emitted 
bosons become incoherent with their source and among each other. Therefore, the coherence 
between the source and the emitted bosons is lost. A complementary argumentation with 
the normalized real-space correlation functions is deferred to the Supplementary Information 
and Fig. S2 therein. 

In the spirit of the seminal work of Penrose and Onsager on reduced density matrices 
[22], we tackle the question "How strong is the loss of coherence in many-body systems?". 
The natural occupation numbers, pf^{t), obtained by diagonalizing the reduced one-body 
density (see Methods for details) define how much the system can be described by a single, 
two or more quantum mechanical one-particle states. The system is condensed and coherent 
when only one natural occupation is macroscopic and it is fragmented when several of the 
pf^it) are macroscopic. In the right panel of Fig. |4] we plot the evolution of the natural 
occupation numbers for the studied systems as a function of propagation time. The initial 
system is totally coherent - all the bosons reside in one natural orbital. However, when 
some fraction of the bosons is emitted, a second natural orbital gradually becomes occupied. 
The decaying systems loose their coherence and become two-fold fragmented. For longer 
propagation times more natural orbitals start to be populated, indicating that the decaying 
systems become even more fragmented, i.e., less coherent. In the few-boson cases, N = 2 
and N = 4, one can observe several stages of the development of fragmentation ^311 [32] - 
beginning with a single condensate and evolving towards the limiting form of entangled 
fragments in the end which means a full "fermionization" of the emitted particles. In this 
fermionization-like case each particle will propagate with its own momentum. For larger A^ 
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the number of fragments also increases with tunnehng time. The details of the evolution 
depend on the strength of the interparticle interaction and number of particles. By resolving 
the peak structure in the momentum spectra of the tunneling systems at different times we 
can directly detect and quantify the evolution of coherence, correlations and fragmentation. 

Finally, we tackle the intricate question whether the bosons are emitted one-by-one or 
several-at-a-time? By comparing the momentum spectra p{k, t) depicted in Fig.[3]at different 
times it becomes evident that the respective peaks appear in the spectra sequentially with 
time, starting from the most energetical one. If multi-boson (two-or-more-boson) tunneling 
processes would participate in the dynamics, they would give spectral features with higher 
momenta which are not observed in the computed spectra depicted in Fig. |3| A detailed 
discussion and a model are given in the Supplementary Information. This model suggests 
that the bosons tunnel out one-by-one. However, the fact that the peaks' heights and 
positions evolve with time, indicates that the individual tunneling processes interfere, i.e., 
they are not independent. The origin behind this interference is the interaction between the 
bosons. We conclude that the overall decay by tunneling process is of a many-body nature 
and is formed by the interference of different single-particle tunneling processes taking place 
simultaneously. 

We arrive at the following physical picture of the tunneling to open space of an interacting, 
initially-coherent bosonic cloud. The emission from the bosonic source is a continuous, 
polychromatic many-body process accompanied by a loss of coherence, i.e., fragmentation. 
The dynamics can be considered as a superposition of individual single-particle tunneling 
processes of a source systems with different particle numbers. On the one hand ultracold 
weakly interacting bosonic clouds tunneling to open space can serve as an atomic laser, i.e., 
emit bosons coherently, but only for a short time. For longer tunneling times the emitted 
particles become incoherent. They lose their coherence with the source and among each 
other. On the other hand, we have shown the usage of tunneling ultracold atoms to study 
the dynamics of ionization processes. Each peak in the momentum spectrum is associated 
with the single particle decay of a bosonic source made of A^, A^ — 1, A^ — 2, etc. particles - in 
close analogy to sequential single ionization processes. These A^ discrete momenta comprise 
a total spectrum - in close analogy to total ionization spectra. 

As an experimental protocol for a straightforward detection of the kinetic energy of 
the emitted particles one can use the presently available single atom detection techniques 
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on atom chips [5] or the idea of mass-spectrometry (see discussion in the Supplementary 
Information). Summarizing, in many-particle systems decaying by tunneling to open space 
the correlation dynamics between the source and emitted parts lead to clearly observable 
spectral features which are of great physical relevance. 

Methods 

Hamiltonian and Units 

The one-dimensional A^-boson Hamiltonian reads as follows: 

H{xi, • • • > a;Ar) = ^ f --— ^ + V{xj) j + ^ Xo6{xj - Xk). 
j=i ^ ' j<k 

Here Aq > is the repulsive inter-particle interaction strength proportional to the s-wave 

scattering length of the bosons and Xi is the coordinate of the i-th boson. Throughout 

this work A = Xo{N — 1) = 0.3 for all considered is used. For convenience we work with 

the dimensionless quantities defined by dividing the dimensional Hamiltonian by , where 

h = 1.05457 ■ io-34HLkg is Planck's constant, m is the mass of a boson and L is a chosen 

sec ' 

length scale. At t < trap is parabolic V{x,t < 0) = ^x^, the analytic form of V{x,t) after 
the opening is given in ref. [31]. 

In this work we consider ^^Rb atoms for which m = 1.44316 ■ 10^^^ kg and = 90.4ao 
without tuning by a Feshbach resonance, where oq = 0.0529 • 10^^ m is Bohr's radius. 
We emulate a quasi- ID cigar-shaped trap in which the transverse confinement is w± = 
2291.25 Hz, which is amenable to current experimental setups. Following [33], the transverse 
confinement renormalizes the interaction strength. Combining all the above, the length scale 
is given by L = 2m[^°L ~ ' 10"^ ^) the time scale by = 1.37 ■ 10~^ sec. 

The Multiconfigurational Time- Dependent Hartree for Bosons Method 

The time-dependent many-boson wavefunction \E'(t) solving the many-boson Schrodinger 
equation i^^^^ = H'^{t) is obtained by the multiconfigurational time-dependent Hartree 
method for bosons (MCTDHB), see refs. [36l [37]. Applications include new intriguing 
many-boson physics such as the death of attractive soliton trains [38], formation of frag- 
mented many-body states [39] and numerically exact double well dynamics [i0l - li2] . Recent 
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optimizations of the MCTDHB, see e.g. jl3l HI] , allow now for the application of the algo- 
rithm to open systems with very large grids (here 2^^ = 65536 basis functions), a particle 
number of up to = 101 and an arbitrary number of natural orbitals (here up to 8). We 
would like to stress that even nowadays such kind of time-dependent computations are very 
challenging. 

The mean-field wavefunction is obtained by solving the time-dependent Gross-Pitaevskii 
equation, which is contained as a special single-orbital case in the MCTDHB equations of 
motion, see ref. [3^ 1^ . To ensure that the tunneling wavepackets do not reach the box 
borders for all presented propagation times the simulations were done in a box [—5, 7465]. 
In the dimensional units we thus solve a quantum mechanical problem numerically exactly 
in a spatial domain extending over 8.29mm (!). 

Many-Body Analysis of the Wavefunction 

With the many-boson wavefunction '^{t) at hand the various quantities of inter- 
est are computed and utilized to analyze the evolution in time of the Rose sys- 
tem. The reduced one-body density matrix of the system is given by p'^'^\x\x']t) = 
(^(t)|*t(a;')^(a;)|^(t))^ where (x) is the usual bosonic field operator creating a bo- 
son at position x. Diagonalizing p^^\x\x';t) one gets the natural orbitals (eigenfunctions), 
and natural occupation numbers pf^ (eigenvalues) from the expression p^^\x\x';t) = 
J2iLiPf'^(t) {'Pf'~^{^'^t)y 'Pi^'^i^jt)- The latter determine the extent to which the sys- 
tem is condensed (one macroscopic eigenvalue) or fragmented (two or more macroscopic 
eigenvalues) [311 |32l |15]. The diagonal part of the reduced one-body density matrix 
p{x,t) = p^^\x\x';t) is the system's density. The first-order correlation function in co- 
ordinate space g^^\x',x;t) = —£i2MsJ^ quantifies the degree of spatial coherence of the 

y/p{x,t)p{x',t) 

interacting system [28| [30] . The respective quantities in momentum space, such as the 
momentum distribution p{k, t) and the first-order correlation function in momentum space 
g^^\k'\k; t) = —£i2Mt£)=^ are derived from p^^Uxlx'; t) via an application of a Fourier trans- 
form on its eigenfunctions. 

In real space the density-related nonescape probability is given by Pnot,pi'^) = Jin P^^' ^^^^ 
(see [34j for the non-hermitian results). The momentum-density related nonescape probabil- 
ity Pnot p{^) is obtained by least-squares fitting a Gaussian function p'^""*'*(fc, t) = Ae~^^^^^ 
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to p{k,t) in the /c-space domain [— oo,0]. A and B are the fit parameters. We then define 
the momentum-density related nonescape probabihty as 
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Figure 1: Protocol of the tunneling process. Top: A generic density p{x,t < 0) (blue line) 
is prepared as the groundstate of a parabolic trap V{x,t < 0) (dashed black line). The trap is 
transformed to the open shape V{x,t > 0) (black line), which allows the system to tunnel to open 
space. Bottom: Sequential mean-field scheme to model the tunneling processes. The bosons are 
ejected from "IN" to "OUT" subspaces (indicated by the red line). The chemical potential //j 
is converted to kinetic energy Ekin,i- the momenta corresponding to the chemical potentials 
ki = y^2mEj.in^i = y/2mpi] i = N, N — 1, l-j^^pear in the momentum distribution, see Fig. [2| 
All quantities shown are dimensionless. 




Figure 2: Many-body tunneling to open space is a fundamentally exponential decay 
process. To confirm that the fraction of atoms remaining in the trap decays exponentially with 
time, we depict the density related nonescape probabilities Pnot,p{'^) ™ '^^S'l ^^'^ Pnot.pi^) ™0" 
mentum space, indicated by the respective solid green and red lines. All quantities shown are 
dimensionless. 
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Figure 3: The peak structures in the momentum distributions characterize the physics 
of many-body tunneling to open space. The total momentum distributions p{k, t) for N = 
101 (top center) and their peak structures for = 2, = 4, = 101, and the respective 
Gross-Pitaevskii solutions, at times ti < ^2 < ^3 < ^4- The broad Gaussian-shaped backgrounds 
correspond to the bosons remaining in the trap, the sharp peaks with positive momenta can be 
associated with the emitted bosons. For = 2 we find two peaks in panel i), for A^ = 4 we find 
three peaks and an emerging fourth peak at longer times, in panel ii). In panel iv) we find three 
washed out peaks for A^ = 101. The corresponding GP dynamics reveals only a single peak for all 
times in iii). The arrows in the plots mark the momenta obtained from the model consideration. 
All quantities shown are dimensionless. 
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Figure 4: Monitoring the coherence of the system. Left: The first order correlation functions 
in momentum space \g'^'^\k'\k]t)\'^ for N = 101 are plotted at t = 0,400,600,700. At t = the 
system is totally coherent, i.e., 1(7^^^ P = 1. At times t > 0, the system remains coherent everywhere 
in A;-space apart from the region around /c = 1, where we find peaks in the momentum distributions. 
The loss of coherence, |(7^^^P ~ only in these regions allows us to conclude that the source 
(trapped) bosons remain coherent at all times while the emitted ones are incoherent. Right: The 
time evolution of the first few natural occupation numbers pf^{t) for = 2, = 4, and A^ = 101 
bosons. The coherence in the systems is gradually lost with time. The systems fragment because 
more and more natural orbitals become populated. All quantities shown are dimensionless. 
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Supplementary Information 



The Many-Body Physics of Tunnehng to Open Space 

Let us commence this Supplementary Information by a model of the physical processes 
constituting the many-boson tunneling of an initially-coherent bosonic systems to open 
space. The processes are schematically depicted in Fig. 1 of the main manuscript. 

Static picture: Basic processes assembling the many-body physics 

The initial and final physical situations in the "IN" and "OUT" subspaces are intuitively 
clear. The totally condensed initial state lives in the "IN" region and is confined by a 
harmonic potential. Therefore, it can be described by a harmonic oscillator-like product 
state, cf. e.g. ref. [SI]. In the final state all the bosons have tunneled out and live entirely 
in the semi-infinite "OUT" region. According to Gaudin, M (1971) Phys. Rev. A 4:386- 
394, the static many-body solution of the one-dimensional bosonic system with short range 
repulsive interaction on a semi-infinite axis can be constructed as a linear combination of 
many correlated (incoherent) states. This imphes that the dynamical final state of our 
system is incoherent. 

To model the steps translating the fully coherent systems to complete incoherence, let 
us first consider the situation in which exactly one boson has tunneled through the barrier 
from "IN" to "OUT" and has no more connection with the interior. The "IN" -system now 
has N — 1 particles and the "OUT" -system has 1 particle. By assuming that no excitations 
have been produced in the "IN" -system, the trapped bosons' energy is exactly reduced by 
the chemical potential fii = — E^~^ ^ see bottom part of Fig. 1 of the main text. Here 
i?* is the energy of the trapped harmonic oscillator product state with the distribution of 
i bosons in the "IN" subspace. We assume that the chemical potential does not depend 
on the number of emitted bosons, because in "OUT" V{x) 0. Let us further ignore the 
inter-particle interaction in the exterior system. Energy conservation requires then, that 
the chemical potential jii of the first boson tunneled from "IN" to "OUT" region must 
be converted to kinetic energy. A free particle has the kinetic energy Eq^j, = ^ - we 



thus expect the first emitted boson to have the momentum ki = J2mEQ^j,^ = ^2m/xi. 
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The above considerations imply that the many-body wavefunction can be considered in a 
locahzed basis \IN;OUT). The process of emission of the first boson in this basis reads 
|A'";0) — >■ |iV — 1; 1*^1). Here the ki superscript indicates that the emitted boson occupies 
a state which is very similar to a plane wave with momentum ki in the "OUT" subspace. 
Now we can prescribe the process of the emission of the second boson as — 1; l'^'^) 
I — 2; 1*^^, 1*^2). By neglecting the interactions between the first and second emitted bosons 
we can define the second chemical potential as 112 — Ei]sr{N—l)—EiN(N—2). Thus, a second 
kinetic energy 2 — A*2 gives rise to the momentum peak at k2 — -\/2m//2- Generally, 

the chemical potentials of the systems made oi N — i and N — i — 1, i = 0, ...,A^ — 1, 
particles arc different, so the corresponding peaks should appear at different positions in 
the momentum spectra. We can continue to apply the above scheme until the last boson 
is emitted |1; l^^i • • • i'^n-i^ _^ |q. j^fci . . . j^fcjv^ ^j^g bottom panel of Fig. 1 in the main text 
indicates the chemical potentials for these one-particle mean-field processes by horizontal 
lines and the processes by the vertical arrows. This simplified mean-field picture of the 
tunneling dynamics is analogous to N sequential processes of ionization, where the energetics 
of each independent process (channel) are defined by the chemical potential of the respective 
sources made of N,N — 1,N — 2, etc. particles. 

Connection of the Model to the Numerical Experiment 

Let us first compute the momenta available in the system oi N — 2 bosons with inter- 
particle interaction strength Aq = 0.3, following ref. [S2]. The difference between the total 
energies of the trapped system made oi N = 2 and N = \ bosons provides ki = 1.106. 
The second momentum associated with the emition of the last boson from the parabolic 
trap gives k2 = 1.000. A similar analysis done for the system of iV = 4 bosons with the 
interparticle interaction strength Aq — 0.1 [Xo{N — 1) = 0.3] gives ki — 1.106, k2 — 1.075, 
ks — 1.038 and k^ — 1.000 for the first, second, third and fourth momentum respectively. 
To relate the model and the full many-body results we draw the momenta estimated from 
the respective chemical potentials in Fig. 3 of the the main text by vertical arrows. The 
agreement between the momenta obtained from the model and the respective ones from the 
dynamics is very good, see the arrows and the peaks in the orange framed plots in Figs. 3i) 
and ii). From this figure it is clearly seen that the later in time we look at the momentum 
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distributions p{k, t), the closer the peaks maxima locate to the estimated results. Moreover, 
our model explains why for N — 101 the peaks are wash out. The chemical potentials of 
neighboring systems made of a big particle number (101 and 100) become very close and, as 

a result, the corresponding peaks start to overlap and become blurred. Nevertheless, they 
are always enclosed by the first and last chemical potentials contributing, see the labels ki 
and /cjv in Fig. 3iv). 

The good agreement between our model and full numerical experiments validates the ap- 
plicability of the emerged physical picture to the tunneling to open space. We continue by 
excluding the possibility that the observed peaks in the momentum spectra can be associated 
with excitations inside the initial parabolic trap potential. This can be done straightfor- 
wardly by calculating the chemical potentials associated with the configurations where one 
or several bosons reside in the second, third, etc. excited orbitals of the trapped system. 
It is easy to demonstrate that the bosons emitted from these excited orbitals would have 
higher kinetic energies resulting in the spectral features with higher momenta. Since the 
computed spectra depicted in Fig. 3 of the main text do not reveal such spectral features 
we conclude that the excitations inside the initial parabolic trap potential do not contribute 
to the tunneling process in a visible manner. 

The above analysis suggests that the overall many-body tunneling to open space process is 
assembled by the elementary mean-field-like tunneling processes analogous to the ionization 
of the systems made of different particle numbers which are happening simultaneously. We 
also are in the position to deduce now that every elementary contributing process is of 
a single-particle type. Indeed, if it were a two-particle process, the kinetic energy of the 
emitted bosons would have been -^2™ ~ \^m^ ^ • ^'^^ large N one can assume that the 
chemical potentials of the first two processes are almost equal, i.e., fif' ^ p=i 2/ii. The 
momentum associated with a two-particle tunneling process would be k^^^ — ^^Amjii — -\/2^i 
- which is far out of the domain where the peaks occur in the exact solutions. 

Tracing the Coherence in Real Space 

Here we complement our study of the coherence in momentum space given in the 
manuscript by its real space counterpart. To characterize the coherence of the tunneling 
many-boson system in real space we compute the normalized real space first order correla- 
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tion function g^^^^x^Xi, t) at various times t for the system oi N — 101 bosons and depict the 
results in Fig. SI. Prom this figure we see that initially the system is fully coherent, namely 
\g^'^\x\\xi;t = 0)p = 1. For t > \g'-^'>{x[\xi;t)\'^ < 1 only in the "OUT" region, indicating 
that only the emitted bosons quickly lose their coherence. In contrast, the source bosons 
living in the interior around Xi = x'l = remain coherent for all times. This corroborates 
our findings from the first order normalized momentum correlation function g^^\k[\ki;t) an- 
alyzed in the main text: the bosons are ejected incoherently from a source, which preserves 
its initial coherence. 

Direct Detection of the Momentum Spectra 

It remains to line out the possible straightforward experimental verification of the emerged 
physical picture. In typical experiments the bosons are ultracold many-electron atoms in 
a very-well defined electronic state. According to the conjectures put forward above, the 
bosons will tunnel to open space with definite kinetic energy. We propose to detect the 
kinetic energy of the emitted bosons by utilizing the techniques and principles of mass- 
spectrometry as schematically depicted in Fig. S2. One can place an ionization chamber at 
some distance from the trapping potential to ionize the propagating bosonic atom suddenly. 
The respective experimental ionization techniques arc presently available, see e.g. Ref. [S3] 
and references therein. The now charged particle will, by application of a static electric field, 
experience a corresponding driving force and change its trajectory. The trajectory of the 
ionized atom or, alternatively, the trajectory of the ionized electron are completely described 
by the respective driving force, the electronic state of the atom and its initial kinetic energy. 
By using a detector capable to detect the charged atom or a photoelectron multiplier for 
the electrons one can monitor the defiection of the ionized particle from the initial direction 
of propagation. The kinetic energy and, therefore, the momentum of the emitted boson can 
be calculated. By this one can detect in situ the momentum spectra p{k, t) corresponding 
to different tunneling times and study the tunneling to open space as a function of time. 

For the few-particle case it is especially interesting not only to obtain the momentum 
spectra, but also to monitor the time-ordering in which the peaks appear, i.e., to monitor 
the time evolution of the momentum peak densities p{k,t). In such an experiment one can 
see whether the signals corresponding to the different ki, i — 1,...,N, will be detected se- 
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quentially, starting from the largest momentum, or they appear to some degree arbitrarily. 
The latter case is a clear indication that the tunneling is a combination of several single 
particle tunneling processes happening simultenously, as we predict. Additionally, this mea- 
surement would be among the first direct observations of the dynamics of the coherence and 
normalized correlations in ultracold bosonic systems. 

Let us summarize. The deterministic preparation of few particle ultracold systems is 
now possible, see Ref. [S4]. Mass-spectrometry is one of the most well-studied techniques 
and working tools available and even more sophisticated detection schemes have been devel- 
oped on atom chips (see Ref. [S5]). The combination of these facihties makes the detailed 
experimental time-dependent study of the tunneling mechanism feasible at present time. 
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Figure SI: The real-space normalized correlation function of the tunneling to open 
space process: t)p is used to measure the spatial coherence in the decaying system of 

N=101 bosons at various tunneling times. White corresponds to |(7^^^P = 1 and black to I^^^^P = 
The red lines in the top left part separate the "IN" and "OUT" regions. Here white corresponds 
to full coherence and black to complete incoherence. In the "OUT" region the spatial coherence is 
lost with time, i.e., I^^^^P ~ on the off-diagonal \g^^\xi ^ xi\xi; The coherence of the source 
bosons is conserved, because in the "IN" part I^^^^P = 1 for all times. See text for discussion. 



24 



Tunneling 
experiment 



Propagation 








o 




4— » 




N 




'c 




_o 










Deflection: 
Electric Field 



P(k,t) 




Detection 



Figure S2: Proposed experimental realization of the momentum spectroscopy of the 
many-boson system tunneling to open space. At some propagation distance from the ex- 
periment (left panel) the bosons are ionized by, e.g., a laser beam (center left). Subsequently, the 
ions/electrons are deflected by a static electric field and counted by a detector (center right). The 
momentum distribution can be obtained as histogram from different realizations of the few- or 
many-boson tunneling process by detection of the deflected particles (right). 
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